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Abstract 

We consider various methods of flow analysis in heavy ion collisions and compare 
experimental data on corresponding observables to the predictions of our saturation model 
proposed earlier We demonstrate that, due to the nature of the standard flow analysis, 
azimuthal distribution of particles with respect to reaction plane determined from the 
second order harmonics should always be proportional to cos2((/) — ^pi) independent of 
the physical origin of particle correlations (flow or non-flow). The amplitude of this 
distribution is always physical and proportional to V2- Two-particle correlations analysis 
is therefore a more reliable way of extracting the shape of physical azimuthal anisotropy. 
We demonstrate that two-particle correlation functions generated in our minijet model of 
particle production Q are in good agreement with the data reported by PHENIX. We 
discuss the role of non-flow correlations in the cumulant flow analysis and demonstrate 
using a simple example that if the flow is weak, higher order cumulants analysis does not 
significantly reduce the contribution of non-flow correlations to elliptic flow observable V2 
in RHIC data. 
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1 Introduction 



Differential elliptic flow in heavy ion collisions is defined as the coefficient of the second harmonic 
of particle distribution with respect to the reaction plane 

Mpt) = (cos2(0p^ (1) 

where (pp^ is the azimuthal angle of the produced particle with transverse momentum px, 
is the azimuthal angle of the reaction plane and the brackets denote statistical averaging over 
all particles with momentum and over different events. According to emerging RHIC data 
differential elliptic flow is an increasing function of pt for small pt, which stops increasing 
and saturates to a constant at — 2 GeV p, |], If. While the hydrodynamic models 

are successful in describing the low-pr increase of V2{pt) i|, they disagree with the high-p^ 
data. Models incorporating jet quenching on top of hydrodynamic expansion also do not agree 
with the high-pT data very well 0. Describing the quark-gluon plasma (QGP) dynamics by 
covariant transport theory requires either extremely high initial gluon density or very large 
parton-parton scattering cross sections in order to fit the high-p^ elliptic flow data 

In our previous paper on the subject |1|] we proposed a model of non-flow particle correlations 
in the initial stages of heavy ion collisions which described V2{pt) saturation data rather well. 
The model was based on particle production mechanism in the high energy regime when gluon 
and quark distribution functions of the colliding nuclei reach saturation . As was suggested by 
McLerran and Venugopalan [|lOl the dominant gluon production mechanism in the early stages 



of the collision is given by the classical field of the nuclei. The field was found at the lowest 
order in as in Extensive numerical simulations of the full solution have been performed 
in |]I2[ while an analytical ansatz for the corresponding spectrum of the produced gluons was 



written down in [13|. At the high energies achieved by RHIC experiments the classical field 
alone can not account for particle production; thus quantum corrections become important. 
For an ansatz of gluon production cross section in AA including the nonlinear evolution of 
see [|15I. 

The essence of our model proposed in [|l| is the following: to estimate the non-flow contri- 
bution to V2{pt) one has to calculate the single and double inclusive gluon production cross 
sections first in the framework of the simple McLerran- Venugopalan model and then include 
the nonlinear evolution effects [l^ in them. Two produced gluons in the double inclusive cross 
section are of course azimuthally correlated with each other This correlation can contribute 
to V2 after being averaged over all particle pairs, which is proportional to the total particle 
multiplicity squared. The latter was related to the single inclusive gluon production cross sec- 



tion similar to how it was done in |jT6[. Two comments are in order here. First of all double 
gluon production cross section of course can not be given by the classical field and is therefore 
not a classical quantity. Calculating it thus corresponds to the first (order a^) correction to 
McLerran- Venugopalan model. Secondly the correlations in the double inclusive gluon produc- 
tion cross section are not just back-to-back, as one would naively expect. One of the main 
advantages of saturation physics is that one does not have to assume that the momenta of 
the produced gluons are large in order to use small coupling and twist expansions like it is 
done in the collinear factorization approach. Saturation calculations include all twists and the 
coupling is kept small by the large saturation scale Qs |10|, [l^, |18[. When the momenta of the 



produced two gluons are not extremely large the correlations are not only back-to-back since it 
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is not required anymore by transverse momentum conservation as some of the momentum can 
be carried away by other soft particles. A calculation of the lowest order double inclusive cross 
section was done in |T9[ showing not only back-to-back (A0 = vr) but also collinear (A0 = 0) 
correlations. 

Since the exact double inclusive gluon production cross section is not known we constructed 
a simple model of single- and double- gluon production [|l| in the spirit of fc-r-factorization 
approach used in |T^. The model successfully described the saturation of f2(pr) at high pt as 



well as centrality dependence of V2{B). The goal of this paper is to discuss compatibility of our 
model with other observables in the flow analysis. 

There are other types of non-flow particle correlations in McLerran-Venugopalan model 



20| , |21| on top of the ones considered in [||]. Those are due to dependence of classical gluonic 
fields on nuclear overlap geometry and may contribute to V2 pO| , [2l[] . (If taken separately these 
correlations can not reproduce the high-p^ behavior of V2{pt) given by RHIC data.) We will 
argue in Sect. 2 that these correlations are higher order in and are therefore parametrically 
suppressed. 

In Sect. 3 we will demonstrate that even in the case of only non-flow correlations the 
distribution of particles with respect to experimentally determined reaction plane averaged 
over many events is proportional to cos 2(0 — in agreement with STAR data [Q. (</> and 

are azimuthal angles of the particle and reaction plane correspondingly.^]) This distribution 
is due to the fact that the reaction plane is also determined from the second harmonic of the 
multiplicity distribution. Therefore the cos 20 shape of the distribution with respect to reaction 
plane does not reflect any physics. At the same time the amplitude of the distribution is given 
by physical correlations, that is by 2 ^2 Jll . 

In Sect. 4 we show that two-particle correlation functions in our model are consistent with 
the data reported by PHENIX To see that, one has to relax the large rapidity interval 
condition employed for simplicity in , which does not apply to the PHENIX detector and to 
flow analysis in general. 

We conclude in Sect. 5 by discussing higher cumulant flow analysis. Higher cumulants 
analysis was proposed in |^ as a way to reduce the contribution of non-flow effects to V2- 
We consider the case when all of standard (2nd cumulant) V2 is due to non-flow two-particle 
correlations, similar to We then calculate the higher order cumulants in this model and 



extract V2 from them. While parametrically the conclusion of |^ still holds, the numerical 
values of non-flow f 2 extracted from the fourth and higher cumulants for RHIC are not signifi- 
cantly smaller than f 2 from the two-particle correlation function in agreement with the results 
of recent STAR analysis . 

2 Different Types of Non-Flow Correlations 

The two-particle multiplicity distribution is given in our model [||] by 

dN dN dN dNcorr , . 



d'^kidyid'^k2dy2 d'^kidyi d'^k2dy2 d'^kidyid'^k2dy2 



^We would like to point out that the true reaction plane angle used in Eq. ([|) is, in principle, different 
from the reaction plane angle "^ji determined from the flow analysis. 
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with the first term in Eq. (Q) given by a product of two single particle distributions due to 
classical gluon fields (disconnected piece) |]TT], |T2|, ^ and the second term corresponding to the 
correlations estimated in |I| (connected piece). Let us make some parametric estimates in the 
framework of classical particle production where the coupling is small as{Qs) ^ 1 and together 
with atomic number A ^ 1 it forms a resummation parameter a'^A^^^ ~ 1 [|1^]. First we note 
that at the leading order single particle distribution is 



dN 
d?k dy 



a. 



(3) 



with 5*^ = vri?^ the cross sectional area of the nucleus. The gluons produced by classical fields 
have typical transverse momentum of the order kx ^ Qs ^ ^/R and thus do not "know" of 
the nuclear overlap geometry. This way the distribution of these gluon modes in Eq. (|^) is 
azimuthally symmetric. Therefore the leading contribution of the first term in Eq. (^ does 
not contribute to V2- However, the soft momentum tail of the single gluon distribution with 
kx ~ 1/-R does depend on geometry of the overlap [^]. If one assumes that large Qs is sufficient 
to keep Us small even for these soft modes one can perform a calculation of these geometric 
effects similar to how it was done in 



of Eq. (i m 



20, |21| . One would obtain the following corrected version 



dN 
d^k dy 



Si 
a. 



1 + 



(4) 



As Q^g ~ a^A^''^ ~ 1 the correction in Eq. (|J) is of order l/a^. To calculate the contribution of 
this correction to elliptic flow observable one can either use Eq. (^ directly in Eq. (|I]) or use 
the definition of f 2 through two-particle correlations |2^, 



{V2) = ^/{^M^^^M). 

In the end one obtains the contribution to V2 of geometric corrections to classical fields 



(5) 



qeom 



5^ 



In contrast the second term in Eq. (H) is parametrically of order ||T| 



dNr. 



d^ki dyi d'^k2 dy2 



Si. 



This can be understood by combining two classical gluon fields and letting them interact, 
corresponds to squaring Eq. (0), dropping one power of Sf (gluons have to be at the 
impact parameter to interact) and multiplying by (interaction). Using Eqs. (^ and ([ 
contribution of correlations in Eq. (|^) to f 2 can be estimated as [Q] 



(6) 



(7) 



This 
same 
i|) the 



V. 



non-flow 



a. 



To see which contribution to f 2 is dominant we observe that 5*^ ~ A^^'^ ~ 1/a 



since alA^^^ 



Plugging this into Eqs. (H) and 



qeom A i non— flow 

gives f 2 ~ and f 2 



a: 



Therefore f 



non 
•2 



-flow 



1. 

> 
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V2^°^ when the couphng is small indicating that two-particle correlations are parametrically 
more important for f 2 than geometrical dependence of classical fields. In the actual RHIC data 
the saturation scale was estimated to be 2 GeV^ corresponding to agiQs) ~ 0.3 which 
means that geometrical effects of |2^, ^ may still be numerically important, though probably 
only in the soft transverse momentum region of fcr ~ 

Finally we should mention that HBT correlations are also present in the model and should 
in principle be added to the right hand side of Eq. (0). Parametrically HBT correlations are of 
the order 

yd) 



d?ki dyi (Pk2 dy2 a 



. .10 



and are much larger than any other correlations mentioned above being comparable to the 
uncorrelated piece. (The estimate of Eq. @ could be obtained by squaring Eq. (^. To derive 
the HBT correlation coefficient one has to consider interference of the particle production 
amplitudes contributing to the first term in Eq. (0) [^.) However these correlations are 
important only when \ki — k2\<l/ R and could be excluded from experimental fiow analysis by 
imposing corresponding cuts on particle momenta. The contribution of HBT correlations to 
the fiow observables was studied in detail in p7[] . 

We should also point out that non-perturbative non-fiow correlations due to jet fragmenta- 
tion, hadronization and resonance decay are not included in our model and could be somewhat 
important numerically in the data. We assume that their effect is subleading. This assump- 
tion can be illustrated by the following observation. For the elliptic fiow observable f 2 we can 
construct an infrared safe definition using Eq. (^) with px weights. Namely one can define the 
following observable 

, \2 E^,jkikj COS 2{^i -(f) j) 

\^2) - V kk- ^ ^ 

where the sum goes over all the particles involved in a fiow analysis and fcj are their respective 
transverse momenta. This observable is infrared safe, in a sense that it is not sensitive to 
hadronization of partons which involves collinear particle splitting and recombination in the 
final state |2^. This implies that fragmentation functions are not needed in calculation of 



this observable and the perturbative calculation similar to the one carried out in []T| would be 
reliable p8|. The observable in Eq. (p!OD is different from the usual definition of V2 with pt 



weights [^, |2^. For the saturation minijet model of [ij the observable (^2) can be written as 



{u2f = ,2 ^,0 / d^hdyid^k2dy2 ^,^^0] , ^1^2 cos2(0i-02) (11) 
{k±) A^iot J d^kidyid^k2dy2 

with {k±) the averaged transverse momentum of the produced particles and Ntot the total 
multiplicity. The typical transverse momentum in the particle spectrum generated by the 



saturation physics is given by the saturation scale, {k±) Qs |]T0|, ^ |1^, |T^. Therefore it 
would be a reasonable approximation to replace kik2 in the integral of Eq. ( [TT| ) by Q"^, which 
in turn would cancel with {k±)'^ = Ql in the denominator yielding 

^ / C^^fc,C^^Ad|/2 j^T. , COS 2(01 - 02) = (^2)' . (12) 

Ntot J d^k^dy^d^k'?,dy2 
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Therefore the elhptic flow observable V2 calculated in the minijet model of |l| is, to a good 
accuracy, equal to the infrared safe observable M2, and is, therefore, insensitive to parton frag- 
mentation physics. Perturbative calculation of |I| is thus sufficient to give a good estimate for 
f 2 and introduction of fragmentation functions would not significantly change the result. 

3 Azimuthal Correlations: Reaction Plane Analysis 
3.1 Simple Example: Autocorrelations 

Let us consider a simple example of scattering events in which N independent particles are pro- 
duced each having a homogeneous azimuthal probability distribution dn/d(f) = {l/N)dN/d(j) = 
1/271. (We set the normalization to 1 for simplicity.) Let us define the reaction plane in each 
event using all particles according to the standard fiow analysis ^9] . After choosing weights 



to be equal to one the reaction plane angle for t>2 analysis is determined by |^ 

y]^-^ sin 2(hi 

tan2v&« = ^ 13 

Ej=i COS20J 

where out of two roots between and 27r one has to take the one with same signs for cos 2"$^ 
and Z^jLi cos 2(f) j. Here the "reaction plane angle" "ifji of course has nothing to do with the real 
reaction plane since all particles are assumed to be produced independent of geometry and of 
each other. Our goal is to calculate azimuthal distribution of one of the particles with respect 
to reaction plane determined by Eq. (p!3D. Of course the correlation of each of the particles 
with the reaction plane in Eq. (0) is due to our inclusion of this particle in the definition of 
the reaction plane and is not physical. It is, nevertheless, instructive to see what the azimuthal 
shape of these autocorrelations would be. 

The distribution of particles with respect to the reaction plane is defined as 



-rN-T—T- 31— ^han2^'ij-— TV 

1 dWR Jo d(pi 002 • • • d(pN V Y.i=i COS 



sm . 



N 



X 



k=l 



cos2 2^^ 

with the distribution function for independent particles 

dn dn dn dn 



^ ( cos 2^^^ cos 20 J (14) 



>ia(p2 ■ ■ ■ acpN CL<Pl «V>2 U^N 



For simplicity we will put = in what follows. Representing 6- and 6- functions as integrals 
we rewrite Eq. ( |T^ ) as 

dC. dr] 1 d(f)2 ■ ■ ■ d(f)N ^ 



dn 



h d"^ 



R 



I / N Ar COS : 

27r 2m r]-ie Jo (27r)^ ^ 
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After integration over angles the expression becomes 



dn 



^1 dm 



R 



1'„=0 



1 

2^ 



d^ di] 1 
271 2m T] — ie 



g-i5sm20i+ir,cos2^i <^ COS 20i 



N-1 



+ {N-2) 



IT] 



(17) 



where the first term in the curly brackets of Eq. (|T^ corresponds to j = 1 term in the sum of 
Eq. (|T^) , while the second term includes the rest of the sum. To evaluate Eq. (0) in the N ^ oo 
limit let us note that since JoW^^ + ff) < 1 for all real non-zero + v"^ and Jo(0) = 1 the 
integrals in Eq. ([T7|) are dominated by small values of ^ and rj. Expanding the Bessel function 
we write 

N{e + 

4 





N-l 








^ exp 






N-*oo 



The second term in the brackets of Eq. ( [T7| ) becomes 



dCdr]_N 

2tx J-oo 27r 27r "2 



exp 



(19) 



where we put the exponent of Eq. ([T^ to one in our large- approximation. Due to Eq. 
the integral in Eq. (p!7|) is dominated by ^ ~ r/ ~ 1/ -\/iV- Therefore expansion of exponent in 
Eq. (|l3) beyond 0th order gives o(l/A^) subleading corrections. 
The first term in Eq. (pTj) gives 



cos . 



2n 



where we have used 



d^ dr] 



-oo 27r 27rz rj — ie 



exp 



(27r)2 



cos , 



(20) 



= P.V.- + 7ri5(r/). 
r] — ie 7] 

Combining Eqs. (|T9D and (|20|) and inserting non-zero r we end up with 

dn 1 



d^i 



R 



(2vr) 



l + ,/-cos2(0-v]/«: 



(21) 



up to o{l/N) corrections. 

We can draw the following general conclusion from Eq. (pT)). If a particle is included in the 
reaction plane angle definition from the n-th harmonic it is correlated to the reaction plane with 
the shape of the correlations' distribution in the large multiplicity limit given by cos n{(f) — r). 
Note that (auto) correlations leading to cos2{/) distribution of Eq. (^) are non-flow correlations 
by their deflnition. 
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3.2 Correlations with Respect to Reaction Plane 



Now we are going to generalize the simple model of the previous subsection to the case of some 
weak pairwise physical correlations between the particles, similar to the correlations considered 
in [|l|. Namely let us consider the case where the particle number 1, which correlations to the 
reaction plane will be studied, is excluded from reaction plane definition but may have some 
physical correlations with the particles 2, . . . ,N contributing to reaction plane determination. 
We thus remove autocorrelation of particle number 1 with reaction plane in accordance with the 
standard method of flow analysis |^ . However, particle 1 can now be physically correlated to, 
say, particle 2. The latter is included in reaction plane definition and is therefore autocorrelated 
to the reaction plane with cos 20 distribution of Eq. (pT]). Our claim here is that it is this 
cos 2(j) autocorrelations distribution, and not the physical correlations, which will be dominant 
in determining the azimuthal shape of the correlations of particle 1 to the reaction plane. In 
what follows we will keep in mind the cuts commonly employed in flow analysis, which were 
recently used by STAR Q]. In only particles with px < 2 GeV were included in the reaction 
plane definition (corresponding to our particles 2, . . . , N) and correlations of particles with 
Pt > 2 GeV (corresponding to our particle 1) to this reaction plane were studied. 

We start by defining a correlated two-particle distribution given for example by minijet 
correlation discussed in and normalized according to 



drirnrr 1 dN, 



corr 



d(j)i d<i)2 N{N - 1) 
The distribution of N particles now becomes 



(22) 



dn dn dn dn , , dncorr dn dn , , 
= _ l_ (]\f _ X) -— . . . |- . . . (23) 

(i02 ■ ■ ■d(f)j<^ d(f)i d(f)2 d(j)i^ d(j)i d(j)2 dcjy^ dcj)^ 

where the factor of — 1 accounts for all possible pairwise correlations of particle 1 with 
particles 2, . . . ,N. In Eq. (^) we omit possible pairwise correlation of particles 2, . . . ,N with 
each other which are of the same order as the correlations shown in Eq. (PB| ) but do not 
contribute to the azimuthal asymmetry of the distribution of particle 1 with respect to the 
reaction plane, giving only an additive small constant. The correlations of more than one pair 
of particles are neglected in Eq. (p3|) since they are suppressed by higher powers of [|l| . 

To determine the distribution of particle 1 with respect to reaction plane we have to sub- 
stitute distribution of Eq. (E3) into Eq. (1^). Evaluation of the first term in the resulting 



expression could be easily done repeating the steps which led to the first term in Eq. (PT|). One 
thus obtains 



dn 



R 



<I'R = 



dC, drj 1 d(f)2 ■ ■ ■ d(f)N dn. 



(27r)2 + i-oo 271 2TTi ri-ie Jo (27r)^-2 



corr 



X ;^cos20je-'^S-2'^^''2'^^+^''S^=2'°'^2'^^ (24) 

We are interested only in j = 2 term in Eq. ( p^ ) since it comes with cos 202 and can give 
non-trivial correlations with the reaction plane. In all the other terms in the sum of cosines 
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the integration over 02 would give an additive constant small compared to the first term on 
the right hand side of Eq. (P^. Assuming that ducorr / d(f)i d(f)2 is an even function of 0i — 02 
only (which is, probably, true for all perturbative spin-independent pairwise correlations and is 
certainly true for correlations considered in we write 



2n 



dUr 



COS 202 



COS 201 



2n 



drir 



cos 2(01 - 02) 



dfi dfi 1 
cos20i27rt;2(l)t^2(2):7— :7— = cos20i — t;2(l) t;2(2), 

a0i (302 2,11 



(25) 



where we have used the definition of flow observable from Eq. (|^) |2^, ^ and substituted 
dn/d(j) = 1/271. ^2(1) and ^2(2) are elliptic flow variables for the particles 1 and 2 correspond- 
ingly. Let us now consider a specific example of flow analysis when the particle 1 represents all 
particles with fixed value of transverse momentum pt while particle 2 represents particles in 
a broad transverse momentum range which go into reaction plane definition (e.g. all particles 
with Pt < 2 GeV as in STAR analysis 0]). Then ^2(1) would correspond to differential elliptic 
flow V2{pt,B), while ^2(2) would be the averaged elliptic flow V2{B), where B is the impact 
parameter of the collision. With the help of Eq. (^), Eq. ( ^4]) becomes 

dn 



H d^a 



R 



<I'r=0 



Jl^ + cos2<l),^{N-l)v2{pT,B)v2{B) 
[211)^ 271 



d^ drj 1 



-00 27T 27Ti T] — ie 
which can be rewritten in the large- iV approximation as 



N-2 



(26) 



dn 



R 



J_ + C0S2^1^NV2{PT,B)V2{B) 

27rr 27r 



X 



d^ drj 1 
27r 27Ti 7] — ie 



exp 



(27) 



Performing the integrations in Eq. (|2^) similar to how it was done in obtaining Eq. 
reintroducing non-zero ^ r we write 



dn 



d^i 



{27ry 



V7TNv2ipT,B)v2iB) cos2( 



R) 



OD and 



(2^ 



Eq. (^) represents the main point of this Section: the shape of the correlations of particles with 
respect to the reaction plane determined by the standard elliptic flow analysis is proportional 
to cos2(0pj, — \l//j) independent of the physical nature of these correlations. Therefore particle 
distribution with respect to reaction plane carries no information about the azimuthal shape 
of physical correlations. 
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Nevertheless, the coefficient in front of the cosine in Eq. (^) is given by the strength of 
physical correlations. To understand the nature of the coefficient in Eq. ( pS] ) let us calculate 
reaction plane resolution in our model. If we divide all — 1 particles defining the reaction 
plane into two subgroups of (roughly) N/2 particles each and determine the reaction plane 
angles in each subgroup independently (\E'i and ^2)1 the reaction plane resolution would be 
given by 



A = \/2(cos2(^i-^2)) 



(29) 



where the averaging is taken over many same-multiplicity events. To determine the strength 
of the reaction plane correlations in our model we note that the sub-event planes are correlated 
in it only due to pairwise correlations between the particles. One particle from plane is 
correlated to plane '^2 via the distribution in Eq. (pSl) . The same particle is correlated to its 
own plane via autocorrelations of Eq. (|21|). To determine correlations between the two 
planes introduced by this one particle we have to multiply the two distributions and average 
over all angles and momenta of the particle (f). The total correlations would be obtained by 
taking this correlation due to one particle to the N/2 power to account for all N/2 particles 
in the sub-event plane definition. Expanding the resulting expression up to the lowest order in 
correlations we end up with 



(Iti tt 
__«l + _iV..(S)=cos2(*,-*.). 

Using Eqs. ( PP| ) and ( PU] ) we can determine the reaction plane resolution as 



(30) 



A 



V2{B), 



(31) 



which is in agreement with 30]. With the help of Eq. (pT]), Eq. (p8|) can be rewritten as 



dn 



R 



i2ny 



[1+ 2v2ipT,B)Acos2{ 



(32) 



Now we can see that, independent of the nature of correlations contributing to f2, particle 
distribution with respect to reaction plane has a cos2(/) shape with the amplitude given by 
differential elliptic flow observable V2{pt, B) times the event plane resolution. In a recent paper 
by STAR 0] the data on particle distribution with respect to reaction plane was reported, 
which was ffited rather well with the ansatz of Eq. ( p2D Eq. ( ^2] ) tells us that since our 
model describes the data on V2{pt,B) rather well ^ it also agrees with the recent data on 
correlations with respect to reaction plane reported by STAR pi. 



4 Azimuthal Correlations: Correlation Functions 

We saw in the previous section that the reaction plane analysis cannot determine the physical 
shape of azimuthal particle correlations. In this section we consider another method of az- 
imuthal correlations analysis - the correlation function method. We believe that this method 
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correctly determines the shape of azimuthal correlations. The two-particle azimuthal correla- 
tion function is defined by |^ 

ry^AA\ dNreal/dA(f) Nmixed ,oq\ 

dJS. mixed/ d/XCp Mreal 

where dN^eai/ dAcj) is the number of particle pairs observed in the same event with a given 
azimuthal opening angle A0 and dNmixed/ dAcf) is the number of particle pairs selected from two 
different events with the same azimuthal opening angle. Nreai and Nmixed are total numbers of 
pairs in the same and in different events correspondingly. Only events with the same multiplicity 
are selected for the analysis of C(A0). In this section we will calculate the contribution of two- 
particle correlations in our minijet model to the azimuthal correlation function. 

The number of particle pairs produced in the same event with a given azimuthal opening 
angle A0 can be obtained by integrating Eq. (0) 

'^^rf = 2-K [ dkikidk2k2dyidy2 ,^^o, , , (34) 
while the total number of pairs is 

Nreai = [ d'^ ki (i^/c2 di/i dy2 „ , , , , . (35) 
J d^ki ayi a^/e2 ay2 

Particle pair where particles are selected from two different events obviously have no physical 
correlations. Therefore their azimuthal angle distribution is trivial and is given by 



dNmixed 1 



Nmixed dAcf) 2n 
Combining Eqs. (|4|), (|5|), and (g) with Eq. (H) yields 



(36) 



dNcor 



(2vr )^ Jdk.k, dk2 k2 dy, dy2 { ..j^ d^dy. + 

^ d k\ d k2 dyi dy2 ^^2^^ d^k2 dy2 ~^ d'^ki dyi d'^k'z dy2 



C( Arh\ — ' " " " - - ^- \a-Ki ayi d'^k2 dy2 d'^k\ dyi d'^k2 dy2 J ('\7\ 



To calculate the correlation function we, therefore, need to know single- and double- inclusive 
particle distributions. 

We calculated the two-particle multiplicity distribution of Eq. (Q) in our previous work |^ 
assuming that there is a large rapidity interval \yi — y2\ ~ l/cts ^1 between the correlated 
particles. It turns out that this approximation preserves all important ingredients of the model 
as long as we are not concerned with the rapidity dependence of elliptic flow and with the az- 
imuthal angular distributions. At the same time it considerably simplifies calculations allowing 
for a simple interpretation of the final result. However, the large rapidity interval assumption 
does not hold in the actual flow analyses. In order to properly describe data one has to re- 
lax the \yi — y2\ ^ 1 condition. Below we intend to compare our model's predictions to the 



correlation function data reported by PHENIX p2[. The (pseudo) rapidity acceptance of the 
PHENIX detector is limited to the interval of 0.7 units in the central rapidity region and the 
Iz/i — I/2I ^ 1 condition certainly does not hold there. 
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Single particle distribution was written in our model [Q as 



dN 2 as 1 I ^2 dxGA dxGA 

d^kidyi Cp S± ^ J dq'^ d{ki — q^y 



d qi , „ -JT-. r^. (3^ 



Relaxing \yi — y2\ ^ 1 condition to \yi — y2\ ~ 1 yields a more general expression for the 
two-particle distribution 

dNcorr Ncol f d^qi f d'^q2 



d'^kidyid'^k2dy2 tc'^CfS± J Q_l J Q2 



f d^qi f d^q2 .2( ^ u u\ 



X 



dxGA dxGA 



dgl 



dq^ 



^(iv I2' hi,k2,yi - 2/2), 



(39) 



i2 



where A{q^, q^, h,k2,yi- 
kinematics (|yi — 2/2! ~ 1) 



?/2) is the two-to-four particles amplitude in the quasi-multi-Regge 
A was found in an impressive calculation performed in [0. We 
refer the reader to the reference |19[ for an explicit expression for A{q^, q^, ki,k2,yi — y2)- In 
the leading logarithmic approximation — ?/2| ~ l/o^s S> 1 this amplitude reduces to 



A{q^, £2' hi.ki.yi -2/2) 



h.1 h.2 



(40) 



This expression for A was used in our model before [|l| yielding only back-to-back correlations 



between the particles. It was argued in that inclusion of the next-to-leading logarithmic 
corrections is essential for understanding the angular distribution of the correlation function. In 
particular, on top of the back-to-back correlations, they introduce correlations at small angles 
A0 between the particles in a pair 



The unintegrated gluon distribution due to the quasi-classical non-Abelian Weiszacker- 
Williams (WW) field {z) of the nucleus is given by |TT|, 



dxGA{,x^ q^ 
dq^ 



d'ze-'^-i J d26Tr(A^^(0)A^^(^)) 



(2vr)2 
7r(27r)^ 



_ 1 ,2oci 2 



a.z^ 



(41) 



where b is the gluon's impact parameter (which we can trivially integrate over in a cylindrical 
nucleus case considered here) and Qf is a certain scale at which nonlinear nature of the gluon 
field becomes evident. In [|[] we argued that one has to include quantum evolution effects in 
quasi-classical formulae in order to describe some important features of the particle spectrum. 
In particular we took into account the fact that evolved gluon distribution scales as a function 
of Pt/Qs only, where Qs is a saturation scale which emerges from the solution to the nonlinear 
evolution equation . We argued in |jl|] that evolution effects can be mimicked by writing the 
argument of the Glauber exponent of Eq. (|41|) in the form 



In ^ 

5 A 



ln(g./A) 



(42) 
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Figure 1: Two-particle azimuthal angle correlation function given by our model compared to 
PHENIX ^/s = 130 GeV data for different cuts in the transverse momentum and centralities. 
We used the following values of parameters: A = 0.15 GeV, A = 197, = 0.3. 



where ^ = q z and A is the infrared cutoff. Thus, after integration over angles in Eq. (^T]) we 
obtain 



dq"^ TT^ Us Jo ^ agZ 

Formulae (|37|), (|38|) , (|39|) , and (^31) , together with the expression for A given in [|l^] comprise 
all necessary theoretical information for calculation of the correlation function C(A0). In order 
to compare our correlation function to PHENIX data we have to integrate in Eq. (|37|) 
over yi and y2 in the rapidity interval from —0.35 to +0.35 (PHENIX detector acceptance) 
and over ki and k2 in the transverse momentum intervals specified by the data analysis of 



22| . The two-particle distribution of Eq. (^) with A from has a collinear singularity at 
small opening angles A(j) due to the contribution of 1 — > 2 gluon splitting [0. It is possible 
that the corresponding increase of the correlation function would be reduced after inclusion 



of higher order perturbative corrections |31| . The small angle particle correlations also receive 
contributions from HBT correlations, resonance decays and hadronization functions. While 
some of these correlations are excluded by appropriate cuts in the experimental fiow analysis, 
some might still remain. To model the effect of these phenomena on the correlation function in 
our model we have imposed a lower cutoff on the invariant mass of the two particles produced 
in Eq. (|39|), with the value of the cut obtained from comparison to the data ||22|. The best fit 



was obtained by requiring that {ki + ^2)^ = 2 {ki fc2Cosh(A?/) — k^ ■ ^2) > 0.07 GeV^. 

Results of our calculations are presented in Fig. |I|. The values of parameters that were 
used are specified in the caption of Fig. ^ Integrations over q^ and q^ in Eq. (|39| ) were cut 
off in the infrared limit of = and q^ = hj 2A. The value of the saturation scale giving 
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the best fit for 0.3 < j9t < 2.5 GeV data is slightly lower than the saturation scale giving 
the best fit for 0.5 < px < 2.5 GeV data {Qs = 0.8 GeV and Qs = I GeV correspondingly). 
This slight discrepancy is probably due to the fact that non-perturbative (soft) effects become 
more important for lower pt data (0.3 < pt < 2.5 GeV) which could be effectively mimicked 
by reducing the saturation scale. The fall-off of our correlation function at the very small 
angles is an artifact of the applied invariant mass cut. More research is needed to completely 
quantitatively understand the correlation function at small opening angles A0. Overall we 



observe a good agreement of our model with experimental data reported by PHENIX |22 



5 Higher Order Cumulants 



A new method of elliptic flow analysis was recently proposed in p3|. It was suggested to 
measure flow using not only two-particle correlations of Eq. (^) ||2^, ^ but also four-, six- and 



other higher order correlation functions. The flow observable can be extracted from cumulants 
formed out of these correlation functions. For instance the fourth order cumulant for elliptic 
flow is defined as 



23 



C2{4} ^ e 



,2j(<^l+02-<^3- 



,2i(0l+</.2- 



2 (e 



,2i(0l-03) 



and in case of only flow correlations is equal to [53 



C2{4} 



(44) 



(45) 



if flow 



The advantage of the higher order cumulant analysis is in the fact that, as argued in [EH 
is larger than non-flow correlations, the contribution of the latter to V2 extracted from higher 
order cumulants is suppressed by powers of particle multiplicity. For instance, if M particles 
are emitted in a collision it can be easily shown that the contribution of non-flow correlations 
to f 2 from Eq. (|^) scales as 1/ \/M. At the same time a similar analysis shows that non-flow f 2 
extracted from the fourth order cumulant of Eq. ( ^5]) scales as i.e., it is suppressed by 

an extra factor of Corresponding data analysis based on Eqs. (^if) and (|^) has been 

carried out at STAR ||^. 

In saturation particle production models one should use the number of participants Npart ~ 
S^Ql instead of M [^,^. For ^ = 130 GeV at RHIC the number of participants for the 
most central collisions is Npart{B = 0) = 344 [jl6[ which gives a suppression factor ^ 
0.23. This factor increases with decreasing multiplicity and, correspondingly, centrality of the 
collision. We are going to give a somewhat more careful estimate of what exactly the suppression 
factor is in the framework of our model of non-flow correlations. 

In principle to estimate the contribution of saturation physics to the fourth order cumulant 



from Eq. (|^) one has to calculate the four particle inclusive cross section in the quasi-classical 
approximation of |]TU|, |TT|, |T2|, |TB|, [T^. This task seems to be extremely difficult. Constructing 
a model similar to what was done for double inclusive cross section is dangerous since at this 
high order the model inspired by factorization approaches may not be valid at all. What we 
are going to do here is to neglect this four particle correlations and estimate the contribution 
of two-particle correlations of Eq. (^ to the fourth order cumulant in Eq. (|4^) . 
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Substituting the distribution of Eq. (0) into Eq. 
C2{4} = 2 



we obtain 



dNcorr dNcorr ^2i(0i+</.2-03-</'4) 



X 



f dN_ dN_ dN_ dN_ _|_ g T iR. dN_ dNcorr _)_ 3 f dNcorr dNcorr 

J d(j>x d(f)2 d(t>3 d<f>4 J d<f>i d<f>2 d^3d(f)4 •> d(f)id(t>2 d<f>3d<p4 y j , j ^^^^^^ 



( r dN_ dN I r 
\J d0i d(f>2 ~^ J 



dNcor 



(46) 



where the integral sign denotes integration over all azimuthal angles to follow. The numerator is 
the same in both terms in Eq. (0). It gives the prefactor in Eq. (|46|) . However the denominators 
are not exactly the same in the two terms in Eq. (0), and this is reflected by the difference 
in the parentheses of Eq. (^). The denominator of the second term is equal to the number of 
pairs squared, while the denominator of the first term is equal to the number of quadruplets of 
particles. 

Expanding the expression in parentheses of Eq. (HBf) to the lowest non-trivial order in 
dNcorr /d4>id(j)2 yields 



C2{4} 



r dNcorr g2i((jii-03) 

J d4>\d<f>-3 

r dN dN 
■J d<f>i d<f>3 



I 



dNcorr 

d(f>id<f>2 
r dN_ dN 
J d<f>i d<f>2 



-8c2{2}2 



dNcorr 

__d4id4>2_ 
r dN dN 
•' d(t>\ d4>2 



-sv2{2y 



I 



dNcorr 

d(l>id(f>2 
r dN_ dN_ 
•' d0i d4>2 



(47) 



where we follow the notation of ^ 
order cumulant C2{2} (see Eq. 



for elliptic flow variable f2{2} extracted from the second 
. Using Eq. (|46| ) we can deduce from Eq. { ^7\j the following 



expression for the elliptic flow extracted from the fourth order cumulant 

_ 1 

r dNcor- 



V2{4} = V2{2} 



r dN_ dN 

J d<j)i d<f>2 



(48) 



To put a lower bound on the suppression factor in Eq. (|4 
always less or equal to one 



dNcorr 

d<f>id4>2 



r dN_ dN 
J d<f>i d(f>2 



> V2{2Y. 



we first note that since cosine is 



(49) 



According to the ^/s = 130 GeV data reported by STAR the elliptic flow for the most central 
collisions is V2 = 1.87±0.25% which after insertion into Eqs. (49) and (^) gives the suppression 
factor 



V2{2} 



> 0.23, 



(50) 



where the precise agreement with the number cited before is coincidental. Using the full ex- 
pression in Eq. ( P^D instead of Eq. (^S|) does not significantly change the result. 

The suppression factor of Eq. ( ^8]) can be estimated in the framework of our minijet model 
nil. The result is 



^2(4} 

V2{2} 



0.5^0.75, 



(51) 
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where the discrepancy is due to infrared cutoff sensitivity of our modeL An exact calculation 
in the saturation framework of |]TU|, |TT|, |r7|, [T3[ should give a more precise estimate. Both 
numbers in Eqs. ( |50D and ( pTj) are consistent with the STAR data presented in Fig. 13 of 0. 
There the ratio of the elliptic flow extracted from the fourth and the second order cumulants is 
in the range of 0.25 — 0.5 for the most central events, where the discrepancy is due to different 
methods of calculating the fourth order cumulant 0. 

Using the same two-particle correlations model we can estimate the ratio of the elliptic 
flow variable extracted from the 6th-order cumulant to the one extracted from the 2nd-order 
cumulant to be 



^^2(6} 



108 



d4>\d<f>2 

r dN dN' 
J dipi d4>2 



0.5^0.75 



(52) 



where we give a numerical estimate of the ratio in our model for the most central collisions. 

Fig. 13 of [§] shows that the ratio of W2{4}/w2{2} increases approaching unity with decreas- 
ing centrality down to 60% centrality. This is in qualitative agreement with our Eq. (^5]), which 
implies that V2{4:}/v2{2} ~ l/Nl^tt ~ l/(5fQ2)i/4 ^j^^^ ^j^g ^atio increases with decreasing 
centrality. 

The steps shown in Eqs. (^Bj), ( ^7] ) and ( P5| ) can be repeated for differential elliptic flow 
V2{pt) analysis in the framework of the same two-particle correlation model of Eq. (H). Fixing 
the transverse momentum of one of the particles in the cumulant analysis to be pt we obtain 
the following expression for differential elliptic flow extracted from the fourth order cumulant 
V2 {4} [pt) in terms of the differential elliptic flow from the second order cumulant V2 {2} (px) 



V2 {4} {pt) ^ V2 {2} (p2 



/ 



dN^, 



I 



dN dN 



1 + 



dNcorr dN 

d'^PTd(j)2 d4 >:i 
dN dNcrn 
d^PT d(p3d4 



(53) 



where the integral sign still implies integration over all azimuthal angles following it. At high 
transverse momentum pt the ratio in the parenthesis of Eq. (|53D goes to a constant since both 
the numerator and the denominator scale in the same way with pt [0]. Therefore, since in our 
model [0 f 2 {2} (px) goes to a constant at high pt, the differential elliptic flow extracted from 
the fourth order cumulant f 2 {4} (p^) would also approach a constant at high px in qualitative 
agreement with the STAR analysis To estimate the value of this V2 {4} (px) asymptotics 
we employ Eq. (O) to write 



V2 {4} (pt) 
V2 {2} (pt) 



(0.5^0.75) 



l + lng,/A 



0.7^1.0 



(high Pt) 



(54) 



for Qs = 1 GeV and infrared cutoff A = 0.15 GeV. The numbers shown in Eq. (|54D appear to 
be consistent with the STAR results shown in Fig. 15 of |Q . 

Despite the qualitative (and semi-quantitative) agreement of our model with the data of 
we do not attempt to flt the STAR fourth order cumulant data at the moment since for a 
consistent flt one needs to include the contributions of the 4-particle correlations in Eq. (0). 
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